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Abstract 

This paper solves the integral equation which describes the oscillating inhomogeneous string, by 
using a spectral expansion method in terms of Chebyshev polynomials. The result is compared 
with the solution of the corresponding differential equation, obtained by expansion into a set of 
sine-wave functions, with emphasis on the accuracies of the two methods. These accuracies are 
determined by comparison with an iterative method which allows a precision of 1 : 10^^. The 
iterative method is based on a old method by Hartree, but contains innovative spectral expansion 
procedures. TgX 
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I. INTRODUCTION 



The teaching of computational physics courses is now practiced by many universities, and 
excellent text books are available supporting this endeavor |l| , as well as papers describing 
such courses [2]. In particular, the vibrating string provides an excellent topic jsl, since the 
solution of the corresponding differential equation can be achieved by several different ways, 
and useful comparison between the different methods can be provided. 

If the string is inhomogeneous, the separation of variables method becomes more involved 
than for the homogenous case, since the spatial part becomes the solution of a Sturm- 
Liouville (SL) eigenvalue equation and is no longer a simple sine wave. The SL equation is 
usually solved by expansion into a basis set of functions (sine waves for the clamped string) 
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that will lead to a matrix equation for the expansion coefficients. The eigenvalues and 
eigenvectors of this matrix then provide the SL functions, but, the accuracy depends on the 
size of the basis, and correspondingly on the size of the matrix. The accuracy of this method 
can be studied by introducing an entirely different method of solution of the SL equation, 
which is normally not discussed in the existing teaching literature. This method, denoted as 
I EM (for Integral Equation Method), consists in transforming the SL differential equation 
into an equivalent integral equation, and solving the latter by an expansion into Chebyshev 
polynomials j4|, jsl. This method has the advantage that its accuracy can be automatically 
pre-determined by means of an accuracy parameter, the number of mesh points required to 
achieve a particular accuracy is much smaller than for the more conventional finite difference 
methods (by a factor close to 20), and the size of the matrices is kept small by a partition 
technique, thus avoiding the drawbacks of large matrices in conventional integral equation 
solution methods. These advantages are important for the solution of a computationally 
complex problem {g]. It is the purpose of the present paper to explain the I EM method in 
simple terms, and apply it to the solution of the inhomogeneous vibrating string. A method 
to solve the SL iteratively, thus avoiding the introduction of the inaccuracies described 

n 

above, will also be presented. This method was ffist devised by Hartree [7| in the solution of 
atomic physics energy eigenvalues of the Schrodinger equation. It has now been adapted to 
the spectral I EM solution of the equivalent integral equation j^, and since it can achieve 
an accuracy of 1 : 10^^ it does provide the bench mark values against which the previous 
methods for the inhomogeneous string can be compared. 

The method for the solution of the string equation is very close to the solution of the 
important quantum mechanical time independent Schrodinger equation. Since the properties 
of the string are much easier to visualize than the properties of the Schrodinger equation, the 
present discussion of the vibrations of the string also serves as a pedagogical introduction 
to the numerical methods required for quantum mechanics. The numerical calculations are 
done with MATLAB. An excellent introduction into both MATLAB and numerical methods 

r 

can be found in the book by Recktenwald [S] . The MATLAB programs for the calculations 



io|. 



presented here will be available in the "compadre" digital library 

In summary, the main purpose of this paper is to introduce to the teaching community 
the use of spectral expansions, especially for the solution of integral equations, because of its 
elegance, its accuracy, and its computational economy. The method of spectral expansions 
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is not new (since circa 1970) and is described in the excellent book by L. N. Trefethen 

Section 2 presents the differential equations describing the vibrating inhomogeneous string 
and the solution by expansion into a basis of sine-wave functions; Section 3 presents the 
basics of expansions into Chebyshev polynomials, section 4 presents the Sturm-Liouville 
(SL) integral equation that is equivalent to the differential equation, and also presents the 
solutions in terms of the lEM spectral method; in section 5 the iterative solution of the 
SL equation is described, and the accuracy of the previous methods is examined. Section 6 
contains a summary and conclusions. 

II. THE INHOMOGENEOUS VIBRATING STRING 

Consider a stretched string of metal, clamped between two horizontal points Pi and P2. 
The distance between the fixed points is L, the mass per unit length p of the string is not 
a constant, as described below, and the speed of propagation of the waves depends on the 
location along the string. When a disturbance is excited along the string, the particles on the 
string vibrate in the vertical direction with a distribution of frequencies to be determined. 

Denote by y{x, t) the (small) displacement of a point on the string in the vertical direction 
away from the equilibrium position y = 0, for a given horizontal distance x of the point from 
the left end Pi, and at a time t. As can be shown, the wave equation is 

dx^ Tdt^ ^ ' 

where T is the tension along the string. We define a function R{x) which is dimensionless, 
and which describes the variation of p with x according to 

p{x) = po R{x) (2) 

where po is some fixed value of p. Defining a reference speed c according to 

(3) 

the wave equation becomes 

^-^R(x)^ = (4) 
dx^ c^^'dt^ ^' 

According to the solution by means of separation of variables, y{x,t) = ip{x) A(t), one 
obtains the separate equations 

+ AR{x) ^{x) = (5) 



and 

'LA . (6) 

We assume that the constant A is positive. A general solution for Eq. ([6]) is acos(wt) + 
hsm{wt)^ with 

w = cVA (7) 

where A is an eigenvalue of Eq. ([5]) . ^ ^ 

The Eq. (|5]) is a Sturm-Liouville equation [12] with an infinite set of eigenvalues 
A„, n = 1,2,3,... and the corresponding eigenfunctions ipni^) form a complete set, de- 
noted as " sturmians" , in terms of which the general solution can be expanded 

oo 

y{x,t) = ^[a^ cos(w„)f:) + 6„ sin(w„t)] i^nix) (8) 

n=l 

where W7„ = c^/]v^. The objective is to calculate the functions ipn{x) and the respective 
eigenvalues A^ as solutions of Eq. (|5]), with the boundary conditions that y = for x = 
and X = L, 

MO) = ML) = 0, (9) 

and that for t = 

y{x,0) = f{x) and dy/dt\t=o = g{x). (10) 

The constants a„ and hn in Eq. ([8]) are obtained from the initial displacement of the string 
from its equilibrium position f{x) and g{x), in terms of integrals of that displacement over 
the functions ilJni.^). 

f{x) V^„(x) dx; hn = — / g{x) ^„(x) dx. (11) 

Jo 

A. The case of the homogeneous string. 

In the case that the string is homogeneous, the function R{x) = 1 becomes a constant, 
and the Sturmian functions are given by the sine functions, i.e., ipn{x) = 4>n{x), with 

4>n{x) = a/2/L sm{knx), kn = n{7i/L), n = l,2,3... (12) 

and the eigenvalues become A„ = fc^ = [nvr/L]^. Assuming that the initial displacement 
functions f{x) and g{x) of the string are given by 

f{x) = xsm[{n/L)x], g{x) = (13) 




0.2 0.4 0.6 0.8 1 



X 

FIG. 1: Vibrations on the homogeneous string. The symbols * mark the initial displacement of 
the string from its equilibrium position, given by Eq.(|13p. The numbers written next to each curve 
indicate the time, in units of L/c 



and 

L = lm, c = 800 m/s. (14) 

then one can evaluate Eq. ( ITT]) for the coefficients a„ analytically (all the 6„ = 0). One finds 
that all Qn vanish for n odd, with the exception for n = 1, for which 



For n even, the corresponding result for a„ is 

1 1 



r2 / 2 

n = 2,4,... (16) 



7r2 V L 



[1 + ny (l-n)2_ 



With the above results the truncated sum ([8]) 

n max 

?/("™^^)(x,t) = J2 COs{Unt) + bn Sm{Unt)] 0„(x) (17) 

n=l 

can be calculated. The result is displayed in Figs. ([T]) and ([2]) 

For n ^ 1, a„ will approach like (1/n)^, i.e., quite slowly. It is desirable to examine 
how many terms are needed in the numerical sum of Eq. (ITTI) in order to get an accuracy of 
4 significant figures in y. A good guess is that the sum of all terms not included in the sum 

y a„ cos(a;„t) ~ -4— W — / —cos{—tn)dn (18) 

nmax + 1 --nmax + i 
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FIG. 2: Continuation form Fig.([T]) of the time development of the vibrations of the string. 



should be less than t/max x 10~^. The integral in Eq. f|T8|) is smaller than Ximax+i(V'^)^'^^ — 
(^max + 1)^^/2 (since the cos term produces cancellations), and one obtains the estimate 

I ^ a„cos(a;„t)| < 2— W — (nmax+1)"^ (19) 

n max +1 

With nmax = 50 the right hand side of Eq. (fT9l) is ~ 10~^. A numerical evaluation of the 
difference \y^^^\x,0) — /(a;)| is less than 10~^, which confirms that with nmax = 50 the 
accuracy expected for y'^^^\x,t) is better than 1 : 10^. 



B. The inhomogeneous string by means of a Fourier series expansion 

An approximate solution to Eq. (j5]) for i/jn is to expand it in terms of the Fourier sine 
waves given by Eq. f[T^ . since these functions obey the same boundary conditions as the 
ip'nS. The approximation consists in truncating that expansion at an upper limit imax = N, 
and also drop the sub- and -superscript (n) for the time being 

N 

i;(^\x) = J2di'Mx). (20) 

£'=1 

Inserting expansion (!20!) into Eq. ([5]), remembering that (P(f)e{x)/dt'^ = —kj(j)i{x), multiply- 
ing Eq. ([5]) by a particular function 4>e{x), integrating both sides of the equation over dx 
from X = to X = L, and using the orthonormality of the functions 4>e{x), one obtains 

N 

-kj df, + K^ Re/'de = (21) 

£'=1 
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where 



Rfj' 



(t>e{x)R{x)(f)e'{x) dx 



(22) 



are the matrix elements of the function R over the basis functions (f)i. This equation f l2T|) 
can also be written in matrix form, where 



A? 


\ 








(Ri,i 


Rl,2 


Ri,3 ■ ■ 






(dA 








d2 




R2,l 


R2,2 


R2,3 ■ ■ 


■ R2,N 




d2 








ds 


= A 


Rs,! 


-^3,2 


R3,3 ■ ■ 


■ R3,N 




d3 


\ 






\dN/ 




\Rn,i 


Rn,2 


Rn,3 ■ ■ 


■ Rn,n / 




\dN/ 



(23) 



or more succinctly 



k^d = A Rd, 



(24) 



where bold letters indicate matrices, and a vector quantity indicates a (A^ x 1) column. Since 
all the kes are positive, the matrix can be defined as 



^2 



V 



(25) 



and one can transform Eq. into 



where 



and 



M 



fourier 



1t-> i.-l 



k^R k 



Ufi k dn ■ 



(26) 



(27) 



(28) 



While Eq. ( 1241) is a generalized eigenvalue equation, Eq. ( 1261) is a simple eigenvalue equation. 
The vectors Un are the N eigenvectors of the N x N matrix M. fourier, and 1/A„ are the 
eigenvalues. Furthermore, since R is a symmetric matrix, NL fourier is also symmetric. The 
eigenvectors of a symmetric matrix are orthogonal to each other, i.e. (uri)'^ ■ Um = 5n,m- Here 
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T indicates transposition. However the vectors dn are not orthogonal to each other, since 

{dn) ■ dm — {Ufi) k Um- 

In summary, the procedure is as follows 

1. Choose an upper truncation limit of the sum ( l20l) : 

2. Calculate the matrix elements Ri^i' so as to obtain the N x N matrix R 

3. Construct the matrix NLfourier from Eq. f l27|) . and find the eigenvalues (1/A„) and 
eigenvectors Un, n = 1,2, ..N, by using the MATLAB eigenvalue command [V, D] = eig(Wl). 
The output D is a diagonal matrix of the eigenvalues and V is a full matrix whose columns 
are the corresponding eigenvectors so that M * V = V * D. For example = V(:, n) 

4. If $(a;) is the column vector of the basis functions (f)e{x), then iIj{x) can be written 
as (the superscript (A^) is dropped now) 

^|Jnix) = (M„/k-i ■ $(x). (29) 

5. In view of Eq. fl29|) the coefficients a„ =and =< g tpn > can be written as 

an =< f^n > = (^„/k-i- < / $(x) > (30) 
bn =< gipn > = {unfk-^- < g $(x) > (31) 



where (/ $(x)) is the column vector of the integrals (/ (f)e) = f{x)(f)e{x)dx, i = 1,2, ..N. 
6. The final expression for y{x,t) can be obtained by first obtaining the coefficients e„ 



Ti -1 



e„(t) = {ur,Y k 



(/ $(x)) cos(w„t) + {g $(x)) — sin(M;„t) 



(32) 



and then performing the sum 



N 



y{x,t) = ^en(t)V'n 



[X 



(33) 



n=l 



In the above, e is the column vector of all e„'s, and \l/ is the column vector of all V'n's- In 
the present discussion we limit ourselves to calculating the eigenvalues A„. 

Assuming that the mass per unit length changes with distance x from the left end of the 
string as 

R{x) = 1 + 2x^, (34) 



9 



and c, /(x) and g{x) are the same as for the homogeneous string case, 



L = Im, c = 800 m/s, f{x) = x sm[{n/ L)x], g{x) = 0. (35) 

then the integrals ( 122|) for the matrix elements Re^e> can be obtained analytically with the 
result 



Re,,, = 2 * 2 ( - ) (-1) 



Rff = l + 2L 



r ^ 1 



1 



^ i' (36) 
(37) 



3 (27r£)2j ' 

The increase of R with x can be simply visualized with the choice flMj) . More realistic 
situations, such as the distribution of masses on a bridge, can be envisaged for future appli- 
cations. 

The numerical construction of the matrices R and Mfourier is accomplished in the MAT- 
LAB program string _f ourier .m which in turn calls the function inh-Str_M.m, using the 
input values 

L = lm, c = 800 m/s, (38) 

The truncation value of the sum Eq. ( I20l) is set equal to either 30 or 60, and the 
corresponding dimension of the matrices Nlfourier or H is N x N. These values are chosen 
so as to examine the sensitivity of the eigenvalues to the size of the matrix Mfourier- 

The results for the eigenvalues A„ are shown in Fig. (|3]) and the corresponding frequencies 
are shown in Fig. The frequencies for the homogeneous string, i.e., for R{x) = 1, are 
shown by the open circles in Fig. (jlj). Since the inhomogeneous string is more dense 
at large values of x than the homogeneous one, the corresponding eigenfrequencies are 
correspondingly smaller. It is noteworthy that the eigenfrequencies of the inhomogeneous 
string nearly fall on a straight line, which means that the frequencies are nearly equispaced, 
i.e., they nearly follow the same harmonic relationship as the ones for the homogeneous 
string. The physical explanation for this property has not been investigated here, but could 
be connected to the fact that the waves for the high indices have more nodes than for the low 
indices, and hence lead to better averaging in a variational procedure. Near the fundamental 
frequency slight deviations from harmonicity do occur, as illustrated in Fig. ([5]). However, 
small deviations from harmonicity will also be caused by other effects such as the stiffness 
of the string. 
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N = 30 




5 10 15 20 25 30 

eigenvalue index 



FIG. 3: The eigenvalues of the matrix Mfourier, defined in Eq. (l27p . The quantity indicates the 
truncation value of the sum in Eq. (|20p . that expands the string displacement eigenfunction ipn{x) 
into the Fourier functions (peix). The dimension of the matrix Mfourier is N x N. 




5 10 15 20 25 30 

eigen-value index 



FIG. 4: The frequencies in units of radians/sec of the vibration of the inhomogeneous string, 
compared with the frequencies of the corresponding homogeneous string. The higher frequencies 
become inaccurate when the dimension of the matrix Mjourier is too small. 

Figures and (jl]) show that for the truncation value of 30, the eigenvalues become 
unreliable for n > 22. This is a general property of the high-n eigenvalues of a matrix, 
which however can be overcome by using the iterative method described further on. The 
table [I and Fig. ([6]) give a quantitative illustration of the dependence of the eigenvalue on 
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5? 0.2 
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I -0.2 
a 
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o 

I -0.6 

I -0-8 

o 

g -1 

Q> 

- -1.2 

5 10 15 

index n 

FIG. 5: The deviation from harmonicity as a function of the eigenfrequency index, for two different 
inhomogeneities. This deviation is defined in terms of the difference between two neighboring 
frequencies d{n) = [w{n) — w{n — 1)] as {d{n + l)/d{n) — 1} * 100. The inhomogeneity is given by 
R{x) = 1 + Fq with Fq either 2 or 4. 



n 


N = 30 


N = 60 


1 

20 


1.614775590198150e-001 


1.6147755902115e-001 


4.092e-004 


4.0933853097811e-004 



TABLE I: Eigenvalues of the matrix M for two different dimensions N x N 
the truncation value by the comparison of two eigenvalues for the same n of the matrix 

Mj^onrier(30 X 30) with thoSe of Mjo„rier(60 X 60). 




-^F0 = 2 
-+- FO = 4 



III. SPECTRAL EXPANSIONS INTO CHEBYSHEV POLYNOMIALS 



First some basic properties of Chebyshev polynomials will be described, then the Curtis- 
Clenshaw method for expanding functions in terms of these polynomials will be presented, 
with special emphasis on the errors associated with the truncation of the expansion, and 
finally the application to solving integral equations will be presented. 
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10 15 20 
eigen-value index 



FIG. 6: The dependence of the eigenvalues of the matrix Mjourier on the dimension N x N the 
matrix. The y— axis shows the absolute value of the difference between two sets of eigenvalues, one 
for N = 30, the other for N = 60. Some numerical values are given in Table [H 



A. Properties of Chebyshev Polynomials 

Chebyshev Polynomials T^{x) provide a very useful set of basis functions for expansion 
purposes . A short review of the main properties needed for the present application 

is presented below. The variable x is contained in the interval —1 — )■ +1, and is related to 
an angle by x = cos 6. This shows that the x's are projections on the of the tip of 

a radius vector of unit length that describes a semi-circle as 6 goes from to vr. In terms of 
the X— variable the T„'s are given by 

To = l 
Ti = X 
T2 = 2x2-1 

T„+i = 2xT„ - T„_i (39) 

In terms of the 6 variable they are given by 

Tn = cos(ra 6); 0<9 <7r. (40) 

It is clear from Eq. ( I40l) that —1 < Tn{x) < 1, and that the larger the index n, the 
more zeros these polynomials have. The T^s are orthogonal to each other with the weight 
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FIG. 7: Chebyshev Polynomials for indices v = 0, 1,2, and 3. They are calculated from Ty{x) = 
cos{v * 6), for the equispaced angles 0. 



function (1 — x'^) The integral X 

/ + 1 /"TT 
T„(x) T^{x) (1 -x^)"^/^ dx = J cos(n^) cos(m^) dO (41) 

has the value if n 7^ m, and the values 7t/2 if n = m ^ and vr if n = m = 0. A plot of 
Ty{x) for V = 0, 1,2, and 3 is shown in Fig. ([7]), which also illustrates that for equispaced 
values of 9 the corresponding values of x are not equispaced. 

The values of x, denoted as ^j, for which a particular T„ = 0, are also not equispaced. As 
can be seen from Eq. (HOi) the zeros of Tn+i{x) with i = 0,1,2, ..N are given by 

■(2i + l) 



6 = cos 



B. The Expansion Method 



2N + 2 



-TT 



0,1, 2, ..AT. (42) 



Given a function /(r), defined in the interval [a, b], in order to expand it into Chebyshev 
polynomials, the first step is to transform the variable r to a new variable x defined in the 
interval [—1, +1]. This can be achieved by means of the linear transformation 

r = a x + /3, (43) 
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with a = {h — a)/2 and /3 = (6 + a)/2. In terms of the x— variable one obtains the function 
/(x) = /(r), and the desired (truncated) expansion is 

TV 

/W(a;) = ^ a„T„(x). (44) 

n=0 

The conventional method of obtaining the expansion coefficients a„ is to multiply Eq. 



on both sides by Tm{x) / yl — x^, integrate over x from —1 to +1, and use the orthogonality 
condition ( 1411) . A more computer friendly alternative was given by Clenshaw and Curtis [l5|. 
It consists in writing Eq. f l4^ + 1 times for the zeros C,o,^i, ■■■^n, of the first Chebyshev 
polynomial T/v+i not included in the sum (jSj), and thus obtain + 1 linear equations for 
the A^ + 1 coefficients, 

N 

n=0 
N 



n=0 



N 



n=0 



which in matrix notation has the form 



c* 



ai 



(45) 



where C is known as the Discrete Cosine Transform. The points are denoted as "support 
points" of the algorithm since the function / has to be known only at these points. The 
elements of the matrix C are Cj j = Tj(^j), and its columns are orthogonal to each other. 
After column normalization, one obtains an orthogonal matrix, and hence the inverse 
can be easily obtained, without the need to invoke a numerical matrix inversion algorithm. 
The matrix C"! is denoted as CMl in the MATLAB program [C,CMl,z] = CjOMI{N) 
available in Ref. lOj. The row vector z contains the values in descending order i = 
N, N — 1, ..0. Inserting the values of Oj, obtained from Eq. (145|) into Eq. (14^ . one obtains 
the value of the truncated function f^^\x) at any point in the interval [—1, +1], and hence 
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the procedure is an interpolation method 16|, ll|. Other cosine transforms also do exist, for 



example one based on the Fourier series expansion method. The method is computationally 
fast, in view of the advent of the FFT algorithms, however a comparison of the spectral 
method with this method is beyond the scope of the present article. 

How good is approximation f l44p to f{x)l If the function is differentiable p times, then it 
can be shown [IT] that 

(46) 

where c is a constant that depends on the p's derivative of /. If the function / is infinitely 
differentiable, then p = oo, and the error (H6l) decreases with faster than any power of N. 
This is denoted as the supra-algebraic convergence of the approximation of f^^\x) to f{x), 
a property also denoted as "spectral" expansion of f{x) in terms of Chebyshev polynomials 
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16|. 



According to Luke [l^. Theorem 2 in Chapter XI, section 11.7 

- ~ a^+iT^+iix) [1 + 2x a;v+i/a7v+2] (47) 

In practice, 

\f(^\x) - f{x)\ < |a^+i| (48) 

This property enables one to pre-assign an accuracy requirement tol for the expansion (IHI) . 
Either, for a given value of A^, the size of the partition of r within which the function /(r) 
is expanded can be determined, or, for a given size of the partition, the value of A^ can be 
determined, such that the sum of the absolute values of the three last expansion coefficients 
aN-2, c^iv-i and a^- is less than the value of tol. 

An example will now be given that shows that, if the function / is not infinitely differen- 
tiable, then the corresponding Chebyshev expansion converges correspondingly slowly. The 
two functions to be expanded are 

f^{r) = r^/^ sm{r) (49) 

f2{r) = r sin(r) (50) 

in the interval < r < vr. While /2 is infinitely differentiable, all the derivatives of the 
function /i are singular at r = 0. The results for the Chebyshev expansions for the functions 
/i and /2 using the Clenshaw- Curtis method are displayed in Fig. ([8]). 
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• sin(r).*sqrt(r) 
o sinr.*r 



5 10 15 

order v of the Cheb. exp. coeff. 



FIG. 8: The Chebyshev expansion coefficients as a function of tlie index for tlie functions /i 
and /2 defined in Eqs. (|49p and (|50p . Since tlie derivatives of tlie function fi have a singularity at 
the origin, the Chebyshev expansion converges more slowly than that of /2, which has an infinite 
number of non-singular derivatives. 



An expansion into a Fourier series of the function /2(r) = rsin(7rr) for [0 < r < 1] is also 
carried out for comparison with the expansion into Chebyshev polynomials. One finds that 
all Fourier coefficients a„ with n = 1, 2.. defined in Eqs. (??) through (1161) for L = 1, vanish 
for n odd, with the exception for n = 1. For n 3> 1, a„ will approach like — 4y^(l/n)^, 
i.e., quite slowly. The absolute value of this result is shown in Fig.(|9]). By comparison with 
Fig. ([8]) one sees that the Fourier expansion coefficients decrease with the index n much more 
slowly than the Chebyshev expansion coefficients. 



C. Integrals based on spectral expansions. 

Given a function /(r), defined in an interval a < r < 6, it is the purpose of this sub-section 
to numerically obtain a spectral approximation to the indefinite integral of this function 




As is done in Eq. 0431] the function /(r) is transformed from the variable r to the function 
f{x) for the variable x C [—1, +1]. Then the integral ( |5Ti) becomes 

:^{r) =^^^h{x) (52) 
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even fourier index n 



FIG. 9: The Fourier expansion coefficients of the function /(x) = xsin(7r x) in the interval [0, 1] in 
terms of the basis functions -^sin(n7r x). The analytic result, given by Eqs. (fT5]) and ([T6|) with 
L = 1, is shown by the symbols *. For odd values of n 7^ 1 they are zero. The solid line represents 
an approximation to \an\ — 4\/^{l/n)^. 



where 



Irix) 



f{x')dx'. 



(53) 



It is desired to obtain the spectral expansion of the approximation to Il{x) 



N 



in-) 



X) 



(54) 



n=0 



where it is assumed that f{x) has been expanded in a series of Chebyshev polynomials, 
as given by Eq. (jS]). In view of the integral properties of Chebyshev polynomials, the 
coefficients 6„ can be expressed in terms of the expansion coefficients a„ of f{x), 



(b \ 

hi 



\bMj 



Sj 



(a \ 



ai 



\aN j 



(55) 



by means of the matrix Sl y^], without loss of accuracy. For the integral 



Ir{x) 



f{x')dx' 



(56) 



an expression similar to (1551) exists, with the matrix Sl rep 
sions for the matrices Sl and Sr exist in the literature j4|. 



aced by Sr. Numerical expres- 



18| and are also available from 
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Cheh., N = 58 


2.43532116647 


Cheb., N = 59 


2.43532116702 


quad, acc = 10^^^ 


2.43532116417 



TABLE II: Integral (|58p obtained with the Chebyshev method, for benchmark purposes 

Ref. [l^ under the name SL_SR. In particular, by noting that Tn{l) = 1 for all n, an 
approximation to the definite integral 3r(r2) = f{r')dr' is given by 

N 



n=0 



(57) 



with an error comparable to Eq. ( HHj) . of the order of |6Ar+i|. The above form of the definite 
integral f l57|) is denoted below as Gauss-Chebyshev quadrature. The existence of Eq. f l55|) 
makes the expansion into Chebyshev polynomials very suitable for the numerical solution 
of integral equations, as will be seen below. 
As an example, the integrals 



r^^'^ sin(r) dr 
r sinfr) dr 



(58) 
(59) 



are evaluated below by using Eq. f loTj) . For comparison purposes Ji was also evaluated using 
the MATLAB integration function quad{@myfun, 0, vr, acc), where acc denotes the precision 
to within which the quadrature result is given. The results are shown in the last line of 
table H] 

If one uses an expansion of the integrand /(r) = sin(r) * r^/^ into a set of Chebyshev 
polynomials, and uses the integral properties of these polynomials by means of the function 
[SL, SR] = SL_SR{N), then for 60 support points one gets an accuracy of 1 : 10^^, but 
the convergence with A^ is slow, as is also the case for the expansion coefficients of /i. The 
values of 3i for two values of A^ are shown in Table HT] below. If, on the other hand, one instead 
uses for the integrand the analytic function /2(r) = sin(r)*r, then the corresponding integral 
converges with A^ much faster, reaching machine accuracy for A^ = 18. These convergence 
properties are displayed in Fig. (fTOj) . where a comparison of the convergence using Simpson's 
quadrature method is also shown. 
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FIG. 10: Comparison of the convergence properties of the Gauss-Chebyshev and the Simpson 
integration procedures as a function of the number of support points. The labels 1 or 2 denote the 
integrals Ii = Jq sin(r) r^/^ dr or I2 = Jq sin(r) r dr, respectively. 



IV. THE INTEGRAL EQUATION FOR THE INHOMOGENEOUS STRING. 

In the previous discussion the Sturm-Liouville functions '?/'„(r), solutions of Eq. were 
obtained by expanding them into a set of Fourier functions (peir), and obtaining the eigen- 
functions and eigenvalues of the matrix Mfourier- This matrix consisted of overlap integrals 
of the inhomogeneity function R{x) with the basis functions (j)^{x). In the present section 
three major innovations are introduced: a) we transform the differential equation ([5]) into 
an integral equation, since the numerical solution of the latter is more stable than that of 
the former, b) we replace the need to do overlap integrals by the Curtis Clenshaw method, 
Eq. fj45l) . of obtaining the expansion coefficients, and c) the basis functions are the Cheby- 
shev polynomials for which the expansion series converges much faster than for the Fourier 
expansions. 

The integral equation that is equivalent to the differential equation ([5]) is 

1 

^(r) = - g{r,r') R{r') %l){r')dr' (60) 

where the Green's function Q{r,r') is given by 

Q{r,r') = ——F{r)G{r') for r < r' (61) 
Ij 

Q{r,r') = ——F{r')G{r) for r > r' 
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and where 

F{r) = r; G{r) = {L-r). (62) 

Both functions F and G obey the equation (PF/dr'^ = 0, d'^G/dr^ = and they are hnearly 
independent of each other. Because of the separable nature of Q the integral on the right 
hand side of Eq. flBU]) can be written as 

^ 1 



g{r,r') R{r') 4j{r')dr' = --G{r) / F{r')R{r') ^{r')dr' 

L Jo 

- yF{r) [ G{r')R{r') ^{r')dr' (63) 

In view of the fact that F vanishes at r = and G vanishes at r = L, and hence 
/o^^('", "t') R{r') il){r')dr' vanishes for both r = and r = L, the functions satisfy the 
boundary conditions. A proof that ■?/'(r) defined by Eq. (!60l) satisfies Eq. ([5]) can be obtained 
by carrying out the second derivative in r of Eq. fl63|) . 

The numerical solution of Eq. (160|) is accomplished by first changing the variable r, 
contained in the interval [0, L], into the variable x, contained in the interval [—1, +1], which 
results in the transformed functions ip{x), Q{x,x'), and R{x'). Expanding the unknown 
solution ip{x) into Chebyshev polynomials 

N 

ipix) = a„T„(x), (64) 

n=0 

as was done in Eq. ( Hil) . then Eq. ( l60l) leads to a matrix equation in the coefficients a„, as 
will now be shown. The coefficients a„ can be placed into a column vector 

a = [ao,ai, ..,aAr]^, (65) 

where T means transposition. The values of ipiC.i) at the support points which are the 
zeros of T/v+i, can also be expressed as a column vector 

^ = mo), ^(6), -mN)]^, (66) 

— * 

and the relation between a and tp, already given in Eq. (H5l) . is 

a = C~^ if, {f=Cd (67) 
Another important relation concerns the integrals 

/X pi 
^{x') dx' and ^r{x) = j 0(x') dx\ (68) 
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where is a function defined in the interval [—1,1], and the corresponding expansion coef- 
ficients a„ are given by a = (p. If ^l,r{x) is expanded into Chebyshev polynomials 

n=N n=N 

^l{x) = and = ^ /3(^)T„(x) (69) 

71=0 n=0 

then the expansion coefficients /3 can be expressed in terms of the expansion coefficients a 
of by means of the matrices S/, and S^j, described near Eq. f l5^ . 

= Sitt and = SrS. (70) 



The matrices C, C ^, S/, and Sr can either be obtained from Ref. 10| or can be found in 



Ref.j^]. Making use of Eqs. f E7|) and f lTU]) one can write the Chebyshev expansion of the 
right and left hand sides of Eq. ( BUj) as 

-^a = M.IEM a (71) 



where 



M/ijM = ^ * C"^ * M3 * Di? * C. (72) 



In the above the factor 1/2, comes from the transformation of coordinates from r to x, and 
where the term L was cancelled by the {i/L) in Eq. flB51) : Di? is the diagonal matrix that 
contains the values of R{^i) along the main diagonal, and M3 is given by 

yi^ = DG*C* Sl*C-^ * DF + DF *C* Sr*C-^ * DG. (73) 

The first (second) term in Eq. fl73l) represents the first (second) term in Eq. fl63l) . = 
diag{F) and = diag{G) represent the diagonal matrices having the values of F{C,i) and 
G{^i) along the main diagonal, the being the + 1 support points described near Eq. 

The explanation for Eq. (!7T!) is as follows: the matrix 'WLjem in Eq- (!^2|) is applied to the 
column vector a, the C in fl72|) transforms the a into the vector '0, the factor DR together 
with the factor DG in (1731) transforms iiito G ® R ® ip (the symbol means that in 
G®R each element of the vector G is multiplied by the corresponding element of the vector 
R, and a new vector of the same length is produced), the additional factor produces 
the expansion coefficients oi G ® R® ip, the matrix Sl or Sr transforms these expansion 
coefficients to the expansion coefficients of the respective indefinite integrals, etc. 
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FIG. 11: Accuracy of the eigenvalues of M fourier and Mi em for various values of their dimension 
N X N. The value of is indicated in parenthesis in the legend. For the Fourier method, N is 
the number of basis functions (j)£ used to expand the Sturm-Liuoville eigenfunctions, and for the 
IEM, N is the number of Chebyshev polynomials used in the expansion, which is also equal to the 
number of support points in the interval [0, L]. The accuracy of the matrix eigenvalues is obtained 
by comparison with a highly accurate result of 1 part in 10^^ obtained by an iterative method. 



A. Results 



After choosing a certain value for the number Njem + ^ of Chebyshev coefficients a 
numerical value of the {Njem + 1) x {Njem + 1) matrix f l72|) is obtained, from which the 
eigenvalues (1/A„), n = 1, 2, ...Niem + ^ can be calculated. The MATLAB computing times 
for the Fourier method for N fourier = 30 and 60 combined using the analytic expressions for 
the integrals needed to obtain the elements of the matrix R is 0.91s, while the computing 
time for the IEM matrix method for all three Njem = 30, 60, and 90 values combined is 
0.75s. Hence the IEM method is comparable in complexity to the Fourier expansion method, 
provided that the overlap integrals (122|) are known analytically. However, a disadvantage of 
the IEM for this application is that some eigenvalues are spurious. Their occurrence can 
be recognized in that they change with the value of Njem, and do not coincide with the 
eigenvalues of Nl fourier- 

The accuracy of these two matrix methods is illustrated in Fig. f llip . It is based on the 
iterative method described below, used as an accuracy benchmark, since it gives an accuracy 
of 1 : 10^^ for the eigenvalues regardless of the value of the eigenvalue index n. Figure f lTTj) 
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shows that the accuracy of the lEM matrix method is considerably higher than the Fourier 
matrix method for the low values of ra, but it is not as monotonic as the latter. The figure 
also shows that the accuracy of both matrix methods depends sensitively on the dimension 
N of their respective matrices M. 



V. THE ITERATIVE METHOD 

This iterative method was introduced by Hartree j3] in the 1950's in order to calculate 
energy eigenvalues of the Schrodinger equation for atomic systems. The method was adapted 
to the use of the spectral expansion method {lEM) and applied to the energy eigenvalue of 
the very tenuously bound Helium- Helium dimer 8[ . The version described below for finding 
the eigenvalues that multiply the inhomogeneity function R, with appropriate modifications 
is also suitable for finding the eigenfunctions for more general SL equations, such as the 
Schrodinger equation [l^. The method is as follows. 

For a slightly wrong value Ai of A there is a slightly wrong function ipi that obeys the 
equations 

f^+A,R(r)Mr) = 0- (74) 

This function does not satisfy the boundary conditions at both r = and r = L unless it 
has a discontinuity at some point rj, contained in the interval [0,L]. To the left of rj the 
function ipi that vanishes at r = is called Yi{r), and to the right of rj it is called k * Zi{r), 
and vanishes aX r = L. Here ^ is a normalization factor chosen such that Yi{ri) = t Zi{rj). 
Both these functions rigorously obey Eq. (FMl) in their respective intervals and are obtained 
by solving the integral equations 

Yi{r) = F{r) - Ai g{r,r')R{r')Yi{r')dr' , 0<r<rj (75) 
Jo 

and 

Zi(r) = G{r) -Ai I g{r, r')R{r')Zi{r')dr' , rj <r <L. (76) 

Jri 

These integral equations differ from Eq. ( !60|) . due to the presence of a driving term F or 
G. However, since the second derivatives of these functions are zero, their presence does not 
prevent that Yi and Zi obey Eq. ( 17i|) in their respective domains. 
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The iteration from Ai to a value closer to the true A proceeds as follows. One multiplies 
Eq. dZZD 

+ AiR{r) Y,{r) = 0, < r < ri (77) 

with ilj{r) and one multiplies Eq. (|5]) with Yi{r), subtracts one from the other, and integrates 
from r = to r = rj. One finds that J^'iYI'^ - f'Yi)dr' = (F^V - ^'Yi)r, = (A - 
Ai) f^' Yiipdr'. Here a prime denotes the derivative with respect to r. A similar procedure 
applied to Zi in the interval [r/, L] yields — { {Z[ip — ip'Zi)rj = (A — Ai) J^' t Ziipdr'. Adding 
these two results and remembering that tZi = Yi for r = rj, and dividing the result by 
4'{rj)iZi{ri) one obtains 

A - Ai = ' , . (78) 

This result is still exact, but the exact function ip is not known. The iterative approximation 
occurs by replacing in the first integral in the denominator by Yi, and by t Z\ in the second 
integral, and by replacing '^/'(r/) in the denominators of each integral by either Yi(r/) or by 
t Zxiji). The final result is 

Aa = Ai + ' , . (79) 

In the above, A was replaced by A2 as being a better approximation to A than Ai, and the 
normalization factor t has cancelled itself out. The iteration proceeds by replacing Ai in the 
above equations by the new value A2. 

The derivatives in the numerator of Eq. ( 1791) can be obtained without loss of accuracy 
by making use of the derivatives of Eqs. fl75|) and f l76|) 



Y[{r) = F'{r) + ^G'{r) [ F{r')R{r') Yi{r')dr' + ^F'(r) / ' G{r')R{r') Yi{r')dr' (80) 
L Jo L Jj. 



and 



Z[{r) = G'ir) + ^G\r) f F{r')R{r') Zi{r')dr' + ^F'{r) f G{r')R{r') Zi{r')dr' (81) 

Jrj ^ Jr 



'rj 

with the result at r = r/ 



and 



Ylirj) = 1 - ^ r'R{r') Yi{r')dr' (82) 
^ Jo 

Z[{ri) = -1 + ^ ^{L - r')R{r') Zi{r')dr' (83) 

^ Jr, 
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n 


An 


n 


A„ 


1 


1.61477559021e-001 


26 


2.42220326385e-004 


2 


4.06257259855e-002 


27 


2.24611142229e-004 


3 


1.81281029690e-002 


28 


2.088546473136-004 


4 


1.02131986136e-002 


29 


1.94699775697e-004 


5 


6.54130338213e-003 


30 


1.81936592475e-004 



TABLE III: Eigenvalues of Eq.® obtained iteratively with Eqs. (f79]) 

In the present formulation the dimensions of A are and the dimension of F, G, Y and 
Z are i, where i represents a unit of length. As noted above, the derivatives with respect 
to r of the functions Y or Z or ip are not obtained as the difference between two adjoining 
positions, but rather as the known derivatives of F and G, together with integrals over Y 
OT Z OT ip according to Eqs. (IHOj) and (IHTj) . In the I EM formulation these integrals can 
be obtained with the same spectral precision as the calculation of the functions Y or Z or 
ip: jsl, hence there is no loss of accuracy either for the evaluation of Eq. (!79|) . or for the 
calculation of A, which can be set to 1 : 10^^. However, it is important to start the iteration 
with a guessed value of A that lies within the valley of convergence of Eq. (1791) . These initial 
values can be obtained, for example, from the eigenvalues of the matrix Is/lfourier described 
above, or from a method described in Ref. {sj. 



A. Results for the iterative method 

Some of the values for A^ obtained to an accuracy of 1 : 10^^ by means of the iterative 
method described above are listed in Table llllt so as to serve as benchmark results for 
comparisons with future methods. The starting values Ai for each n are the results of the 
Fourier method described above with = 60. The iterations were stopped when the change 
A2 — Ai became less than 10~^^ (usually three iterations were required), and tol = 10^^^. 

The error of the functions Y and Z is given, according to Eq. ( l48l) . by the size of the high 
order Chebyshev expansion parameters. For the tol parameter of 10~^^ their values stay 
below 10^^^, as is shown in Fig. (fT2|) . Since there is no loss of accuracy in evaluating the 
various terms in Eq. ( I79|) . the error in the iterated eigenvalues A is also given by Fig. (fT2|) . 
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FIG. 12: The y-axis shows the absolute value of the mean square average of the three last Chebyshev 
coefficents in the expansions of the functions Y and Z. As discussed in the text, the error of the 
eigenvalues A is also given by the y— axis. The number of expansion Chebyshev polynomials 
was increased adaptively as the eigenvalue index n increased. The "jumps" in the values of these 
errors is due to the transition from one value of A^ to a suddenly larger value, as is explained in 
the text. 



In order to achieve this type of error, the number of Chebyshev polynomials used for 
the spectral expansion of the functions Y and Z for the solution of their respective integral 
equations was increased adaptively by the computer program. It was found that for n = 1, 
N = 16; for n = 2 to 6, A^ = 24; for n = 7 to 23, A^ = 24; and for = 18 to 30, A^ = 54. 
This procedure of increasing A^ is different from the procedure used in Ref. [sl, where A^ was 
kept constant and the number of partitions was increased adaptively. The latter method 
was required because of the long range (3000 units of length) of the He — He wave functions. 



VI. SUMMARY AND CONCLUSIONS 



The main aim of this paper is to introduce the spectral expansion method for solving in- 
tegral equations to the teaching community, with the hope that this method can be included 
in computational physics courses in the future. Such expansions converge rapidly with high 
precision, and complement the usual finite difference methods in common use today. The 
example used for the application of such a method is the analysis of the vibration of an 
inhomogeneous string in the separation of variables formalism. The spatial basis functions 
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i^n{f)-i n — 1,2,..., form a complete Sturm- Liouville set, the calculation of which is per- 
formed by means of three methods. In method 1 the function is expanded into a basis 
set of sine waves, and the eigenfrequencies and expansion coefficients for each ipn are the 
eigenvalues and eigenvectors of a matrix Nlfourier- In method 2 the differential equation for 
ipn is transformed into an integral equation of the Lippmann Schwinger type, the unknown 
function is expanded into Chebyshev polynomials, and the expansion coefficients are again 
the eigenvectors of another matrix M/bm- The comparison between these two methods il- 
lustrates the differences and advantages of each, especially their properties as a function 
of the size of the expansion basis. In method 3, which has not been presented previously, 
the differential equation for the Sturm-Liouvillc cigenfunction is solved itcrativcly, and the 
auxiliary functions required for the iterations arc obtained from the solutions of the corre- 
sponding integral equation. The advantage of method 3 is that the precision of both the 
cigenfunction and the eigenvalue can be predetermined by specifying the value of a tolerance 
parameter, and further, no eigenvalue calculation of big matrices is required. In the present 
application the results of method 3 were accurate to 1 : 10^^. 

The comparison between the accuracies of the various methods is illustrated extensively 
by means of appropriate graphs and tables. Applications of these methods to other problems, 
such as the solution of the Schrodinger equation, or the heat propagation equation, or 
diffusion equations in biology, are of course quite possible in spite of the present focus on 
the inhomogeneous string equation. 
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